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^. ! Numerical simulation of plasma turbulence in the Large Plasma Device (LAPD) 
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[Gekelman et al, Rev. Sci. Inst., 62, 2875, 1991] is presented. The model, imple- 
mented in the BOUndary Turbulence (BOUT) code [M. Umansky et al., Contrib. 
Plasma Phys. 180, 887 (2009)], includes 3-D collisional fluid equations for plasma 
density, electron parallel momentum, and current continuity, and also includes the 
Ph. effects of ion-neutral collisions. In nonlinear simulations using measured LAPD den- 

O ' sity profiles but assuming constant temperature profile for simplicity, self-consistent 

• i-H ' 

£>y evolution of instabilities and nonlinearly-generated zonal flows results in a saturated 

Q_i! turbulent state. Comparisons of these simulations with measurements in LAPD plas- 



mas reveal good qualitative and reasonable quantitative agreement, in particular in 
frequency spectrum, spatial correlation and amplitude probability distribution func- 
tion of density fluctuations. For comparison with LAPD measurements, the plasma 
density profile in simulations is maintained either by direct azimuthal averaging on 



each time step, or by adding particle source/sink function. The inferred source/sink 
values are consistent with the estimated ionization source and parallel losses in LAPD. 
These simulations lay the groundwork for more a comprehensive effort to test fluid 

H ' turbulence simulation against LAPD data. 
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I. INTRODUCTION 

Turbulent transport of heat, particles and momentum has an impact on a wide variety of 
plasma phenomena^- 3 -, but is of particular importance for laboratory magnetic confinement 
experiments for fusion energy applications 4-7 . A large number of advances in understanding 
plasma turbulence have been made using analytic theory, for example nonlinear mode in- 
teraction^, instability saturation and secondary instabilities^ -^ 1 -, cascades^ 2 - 1 ^ 3 - and the role 
of sheared flow^— . However it is increasingly necessary to use direct numerical simulation 
as a tool to gain understanding into the complex nonlinear problem of plasma turbulence. 
Additionally, numerical simulation is key to the development of a predictive capability for 
turbulent transport in fusion plasmas. An essential aspect of the development of this capa- 
bility is its validation of numerical simulation against experimental measurement^ 7 -* 1 ^. 

While ultimately validation against measurements in high-temperature fusion plasmas in 
toroidal geometry must be undertaken, it is desirable to have a hierarchy of experiments for 
comparison, with the goal of isolating important physical effects in the simplest possible ge- 
ometry 17 . Linear plasma devices such as LAPD 19 , CSDX-^, VINETA 21 , LMD 22 , HELCAT-2 3 -, 
and MIRABELLE 2 - offer an opportunity to validate turbulence and transport simulations 
in simple geometry and with boundary conditions and plasma parameters with reasonable 
relevance to tokamak edge and scrape-off-layer plasmas. Thanks to their low temperature, 
these devices are highly diagnosable, providing for detailed comparison against code pre- 
dictions. As these plasmas tend to be fairly collisional, fluid (including gyrofluid) models 
have been compared in recent studies, for example, on LMD 25 , CSDX 2 -, and VINETA—"—. 
These studies have not simply compared code output to data, but more importantly have 
been used to extract physics understanding: the importance of ion-neutral collisions in zonal 
flow damping was explored in LMD 25 ; simulations of the VINETA device were focused on 
exploring the formation and propagation of turbulent structures^ 7 -™— ; and recent simulations 
of the LAPD plasma suggest that sheath boundary conditions in some regimes could drive 
strong potential gradients and in this case the Kelvin-Helmholtz instability can dominate 
over drift-type instabilities 30 . 

This paper presents modeling of turbulence and transport in the Large Plasma Device 
(LAPD) 1 ^ using a Braginskii fluid model implemented in the BOUndary Turbulence (BOUT) 
code 3 ^^. LAPD provides a unique platform for studying turbulence and transport. Large 



size perpendicular to the magnetic field (100 < pi/a < 300) results in a large number of 
linearly unstable modes^ and broadband, fully developed turbulence is observed^. Due to 
its length (17m), perpendicular transport can dominate over parallel losses and changes in 
turbulent transport can have a strong impact on radial plasma profiles^. The LAPD plasma 
is similar to tokamak SOL in the sense that the radial plasma density and temperature pro- 
files are determined by the competition of the radial turbulent transport, parallel streaming, 
and volumetric sources. The use of the BOUT code also provides a unique opportunity to 
directly test in linear geometry the same code that is also used to simulate tokamak edge 
plasmas 

Numerical simulations reported here are done in LAPD geometry using experimentally 
measured density profiles. In order to simplify these initial studies, a flat temperature 
profile is assumed, flow profiles are allowed to freely evolve in the simulation, and periodic 
axial boundary conditions are employed. The simulations show a self-consistent evolution of 
turbulence and self-generated electric field and zonal flows. The density source/sink required 
to maintain the average density profile close to the experimental profile is consistent with the 
ionization source and parallel streaming losses in LAPD. Overall, these calculations appear to 
give a good qualitative and reasonable quantitative match to experimental temporal spectra 
and are also consistent with the measured spatial structure and distribution of fluctuation 
amplitude. These results lay the foundation for proceeding with the more difficult task of 
simulations with matched density, temperature and flow profiles along with more realistic 
axial boundary conditions. 

This paper is organized as follows. In Section [TT] the main parameters of LAPD are de- 
scribed, as well as the fluid equations implemented in the BOUT code that are used to model 
the LAPD device. Section II III discusses the two methods of average profile control that are 
used to maintain the average density close to the experimental values. A detailed com- 
parison of simulated turbulence characteristics to the experimentally measured quantities is 
presented in Section [IVl Section [V] discusses the particle transport and the density source 
in the simulations, and also briefly introduces the numerical diagnostics used for verifying 
the solution. Conclusions are presented in Section |VTJ The appendices demonstrate the 
derivation of the azimuthal momentum equation (Appendix |A]) , derivation and discussion 
of the ion viscosity term (Appendix [B]) and a numerical scheme used to avoid unphysical 
solutions due to parallel discretization (Appendix |C]) . 



II. PHYSICS MODEL 

The LAPD device is a long cylindrical plasma configuration with length L ~ 17 m, 
vacuum vessel radius r s =50 cm, typical plasma radius a ~ 30 cm, electron density n e0 < 
5 x 10 12 cm -3 , electron temperature T e < 10 eV, and ion temperature Tj < 1 eV; with an 
externally imposed axial magnetic field magnetic field B z < 0.25 T. Plasmas are typically 
composed of singly ionized helium although argon, neon and hydrogen plasmas can also be 
studied. 

For the calculations discussed here, LAPD is modeled as a cylindrical annulus with inner 
radius 15 cm and outer radius 45 cm. Using the annulus topology allows the LAPD geometry 
to be described in the BOUT code without any modification of the code itself through using 
the built-in tokamak geometry but changing the metric coefficient values^ 3 -. In this setup, 
the poloidal magnetic field of the tokamak configuration corresponds to the axial field of 
LAPD, and the toroidal field is set to zero as it corresponds to the azimuthal direction in 
LAPD. The annulus configuration also avoids the potential complications of the cylindrical 
axis singularity. The magnetic field is taken uniform, along the cylinder axis, and the axial 
boundary conditions are taken to be periodic. 

The simulations presented here are based on the Braginskii two-fluid model 36 . As dis- 
cussed in a linear verification study 33 that uses the same model, collisions play a very impor- 
tant role in LAPD plasmas. Electron-ion collision rate is much higher than the characteristic 
drift frequencies, v ei ^> w*; and the electron mean free path is much shorter than the parallel 
length of the device, A e j <C L\\. Therefore for low frequency, long parallel wavelength modes 
a collisional fluid model might be reasonable choice for modeling LAPD plasmas. Kinetic 
effects can, however, be very important in LAPD, in particular for Alfven waves, where the 
electron thermal speed is comparable to the phase speed of the wave, v^ ~ t>th, e — • Because of 
the large parallel size of LAPD (and plasma beta of order the mass ratio, /3 ~ m e /M), drift 
waves couple to Alfven waves— and kinetic effects may be important 39 . It can be argued 
that even in this case strong collisions may disrupt kinetic processes such as Landau damp- 
ing and fluid description may still provide a good description of the plasma^. Nonlinear 
BOUT simulations using fluid equations can help to identify the limits of the validity of the 
collisional fluid model in LAPD. 



For the simulations described here the following set of equations are used: 

d t N = -v E • ViV - V,| (v l[e N) (1) 

d t v\\ e = -Ve ■ Vv|| e - fi-^V\\N + /iV||0 - u e v\\ e (2) 

8 t w = -v E • VtJ7 - V||(M;|| e ) + b x VN • Vv £ 2 / 2 - Vintx + fkVjZU (3) 

Here A" is the plasma density, v\\ e is the electron fluid parallel velocity, and w is the potential 
vorticity introduced as 

w d = f v ± • {NV ± <f)) (4) 

elsewhere^ except the viscosity term that is added for nonlinear calculations and is discussed 
in Appendix [B] All the quantities here are normalized using the Bohm convention. The 
model used in BOUT is similar to that employed in other efforts to simulate linear devices, 
in particular on LMD 25 , CSDX 26 , and the recent work by Rogers and Ricci in simulating 
turbulence in LAPE>2£. 

Density, temperature and magnetic field are normalized to reference values n x , T ex (the 
maximum of the corresponding equilibrium profiles), and B , the axial magnetic field. Fre- 
quencies and time derivatives are normalized to Q C { X = eBo/m,iC: dt = dt/^ldx, = u/Q C i X ; 



velocities are normalized to the ion sound speed C sx = y/T ex /m,i; lengths - to the ion sound 
gyroradius p sx = C sx /Q C i X ; electrostatic potential to the reference electron temperature: 
= e<fi/T ex . In Eqs. dHHJ) and further, the " " " symbol for dimensionless quantities is 
dropped for brevity of notation. 

While the variables N, v» e and w are advanced in time, equation (JIJ) is solved to recon- 
struct the perturbed potential <fi from w. In the code version used in this work, Eq. (J3J) is 
linearized to increase for computational efficiency, since this equation has to be solved for 
each evaluation of the right-hand side of Eqs. (fUEJ). The vorticity evolution equation (J3j) 
replaces the current continuity equation. Derivation and discussion of this form of equation 
is presented elsewhere^. Derivation of the viscosity term is discussed in Appendix [Bl 

Equations dHHJ) do not include perturbations of the magnetic field. While electromagnetic 
effects are essential to capture the physics of Alfven and drift-Alfven waves, preliminary 
linear simulations indicate that the effect of magnetic perturbations on frequencies and 
growth rates for low-frequency drift wave is small for the LAPD parameters considered 
here. A more detailed nonlinear study of the electromagnetic physics is a subject of future 
work and is outside the scope of this study. 



Time evolution equations (JTJSJ) are implemented in the numerical code BOUT— >22. Orig- 
inally developed for simulations of the tokamak edge plasma, the code has been adapted 
to the cylindrical geometry of LAPD. For the present study special care is taken to avoid 
spurious numerical solutions due to discretization in the coordinate parallel to the mag- 
netic field (see Appendix [U]). Prior to the turbulence calculations presented here the code 
has been successfully verified for a range of linear instabilities potentially existing in the 
LAPD plasma, including the resistive drift, Kelvin- Helmholtz and rotational interchange 
instabilities 33 . 

III. TURBULENT TRANSPORT AND AVERAGE DENSITY PROFILE 

A. Average and local fluctuating fields 

In turbulence where the eddy size is much smaller than the macroscopic system size, 
the separation of spatial scales usually leads to separation of time-scales for the evolution 
of local and spatially averaged fields. In spite of this separation of scales, the average and 
fluctuating quantities are certainly coupled since gradients provide the source of free energy 
driving turbulence; on the other hand, turbulent transport, along with sources, leads to 
evolution of the macroscopic profiles. If no sources are present in the simulation, the profiles 
relax to smaller gradient as is shown in Fig. [TJ For comparison of the simulated and measured 
turbulence characteristics, this profile evolution is usually undesired, since this comparison 
requires collecting a large statistical sample of data for stationary, experimentally relevant 
profiles. 

For the purposes of considering algorithms for average profile control in BOUT, it is 
convenient to represent the fluctuating variables at any spatial location as a sum of the 
time-average and perturbation, 

/(x,£) = ;(x) + /(x,£), (5) 

For cylindrically symmetric configuration it is also useful to separate / into azimuthally 
symmetric and asymmetric components, 

/(x,t) = (/(x,t)) + {/(x,t)}, (6) 

where (/) = (l/27r) f fdO and {/} = / — (/) is the residual. Based on the ergodic 
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FIG. 1. (Color Online) Relaxation of the density profile in a simulation without sources. 

hypothesis, it is assumed that the time-average / is equal to the azimuthal-average (/) , and 
statistical moments of/ and {/} are equal. This separation of variables into an axisymmetric 
and non-axisymmetric part does not preclude a full nonlinear solution in BOUT, but it allows 
easier control over the average profiles of density, temperature and other quantities. 

B. Profile maintenance: suppressing the azimuthal average 



Following self-consistent time evolution of turbulence and macroscopic transport may 
be difficult because the time-scale separation can make such calculations too large to be 
practical. Additionally, including first-principles-based source terms, e.g., for density and 
temperature, can be complicated, involving models for the plasma source, atomic physics, 
radiation transport, etc. 

Without attempting a self-consistent time evolution of turbulence and macroscopic trans- 
port one can consider intermediate time-scales t <^ t ^ T, where the macroscopic profiles 
can be taken as "frozen" , based on known measured experimental average profiles. In this 
case the time-evolution of only the non-axisymmetric part is considered, and a simple tech- 
nique of maintaining the desired average profile is filtering out the axisymmetric part of 
evolving fields. This is illustrated in Fig. ([2}(3]) showing the general appearance of 0, Srii, 



and the evolution of the density and potential fluctuation RMS in a simulation with "frozen" 
density profile. In Fig. (J3J), the potential is split into the turbulence-generated axisymmetric 
component (0), and the non-axisymmetric residual {0} = <p — (0). One can observe the 
development of a zonal flow component (<fi(r)) corresponding to sheared azimuthal flow. 
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FIG. 2. (Color Online) Spatial structure of the perturbed density (left) and turbulence-generated 
potential (right) at t = 5.2 ms. 

The easiest way to control the average profile is to suppress the evolution of the axisym- 
metric component by subtracting the azimuthal average of the right-hand side of Eqs. ([UO]); 
for example, for the density: 

d t N = RHS-{RHS} e (7) 

This is effectively introducing a time-dependent source/sink function necessary to maintain 
exactly the target density profile. 

However, suppressing the axisymmetric part of fluctuations may interfere too much with 
the solution, and one can consider doing this not in the full domain but only on the boundary. 
This would constrain the boundary values, and may be enough to maintain the average 
profile close to the desired. 

C. Profile maintenance through adding sources 

A more physical method to control the profile is to use a source/sink term S(r) designed 
in such a way that the average profile is maintained close to the experimentally measured 
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FIG. 3. (Color Online) Time evolution of the density and potential fluctuation RMS in a typical 
simulation. The potential is split into the axisymmetric ((</>), zonal flow component) and the 
non-axisymmentric (</> — ((/>)) components. 



profile. Rather than developing this source from first principles, an ad hoc source/sink is 
chosen in order to achieve the desired steady state profile in the simulation. As a first 
step, a simulation is performed using the method of subtracting out the azimuthal average 
to maintain the (iVj) profile close to the "target" density profile N i0 (r), which is based on 
a representative experimental probe measurements in LAPD. Once a steady-state turbu- 
lence solution is obtained, the radial particle flux is calculated from fluctuating density and 
potential, 



T = (NiVEr) 



(8) 



Next, the effective volumetric density source S(r) is calculated as 



s = v-r 

9 



(9) 



that can now be added to the density evolution equation, Eq. ([TJ): 

d t N = RHS + S(r) 



(10) 



A subsequent simulation is run with this new source term and without suppressing the 
azimuthal average, allowing turbulent transport to compete with the source term. If needed, 
another iteration can be made by adjusting the source term to account for the mismatch 
between the BOUT predicted profile and the target. A more comprehensive approach to 
self-consistent time-evolution of turbulence and average profiles can be based on adding an 
adaptive source^; however it is beyond the scope of this paper. In the present study extra 
iterations were not necessary; a single step was enough to produce stationary turbulence 
with average density profile close to the target profile, as shown in Fig. HI The evolution of 
the density and potential in a typical BOUT simulation with density source and fixed values 
of the density at the radial boundary is presented in Fig. [5] (animation online). 
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FIG. 4. (Color Online) Instantaneous density profile (JVj) in simulation with density source S(r). 

IV. COMPARISON WITH LAPD DATA 

Before any attempt is made to construct a full simulation of LAPD that self-consistently 
incorporates transport, first principles sources/sinks and profile evolution, it is necessary to 
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FIG. 5. (Color Online) Evolution of density and potential in BOUT simulation with density source. 
Left: density fluctuations; right: RMS of the density perturbations and the axisymmetric and non- 
axisymmetric components of self-generated potential (top); average radial profiles of density and 
potential (bottom) (enhanced online).[URL: |http://dx.doi.org/10.1063/1.3527987.1| 



ensure that the basic characteristics of the turbulence are correctly captured by the phys- 
ical model being used. Initial simulations have been done using periodic axial boundary 
conditions and using the experimentally measured density profile. A constant temperature 
profile (5 eV) is used and the potential is allowed to evolve self-consistently (potential and 
flow profiles are not matched to experimentally measured profiles). The main experimen- 
tal dataset used in this work is taken from Carter et al..—, and includes measurements of 
plasma and turbulence profiles and two-dimensional correlation functions. Density fluctua- 
tions in BOUT are compared against ion saturation current Ii >sa t fluctuations measured in 
the experiment (making the currently unverified assumption that temperature fluctuations 
are negligible). 

To be able to compare turbulence characteristics to the measured data for the experi- 
mentally relevant plasma parameters, the average profiles have to be maintained close to 
the experimental values during the simulation. To do so, the two methods described in 
section IIHI are applied, either subtracting the azimuthal average of the right-hand side of 
the density equation, or adding a time-independent source function S(r) to the right-hand 
side of Eq. (CQ). 

Using these two methods, a steady-state turbulence for LAPD parameters is simulated, 
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solving Eqs. ([UH]) with the average density close to the experimental profile, for a range of 
ion-neutral collisionalities ivinj^ci = 2 x 10~ 4 , 1 x 10 _3 ,2 x 10 -3 ). The estimated value for 
LAPD, based on neutral density n n ~ 5 x 10 n cm -3 , is v in /VL ci = 2 x 10~ 3 . The simulations 
are performed in the radial interval 0.15 < r < 0.45 m, in an azimuthal segment of a 
cylinder of 7r/4 angle, assuming periodicity in the azimuthal angle and parallel direction. 
The grid size is 50 x 32 x 32 points for radial, azimuthal and parallel coordinates. In order 
to improve the statistics, a series of uncorrelated runs is made with slightly different initial 
perturbations. In the experiment, a slow evolution of the mean density is observed on a 
~ 1 — 2 ms time scale. BOUT simulations with fixed background profile can not capture 
the slow variation of the average density because the azumuthally symmetric component of 
density perturbation is continuously removed in BOUT simulation to maintain a constant 
profile. For the purposes of direct comparisons between BOUT results and the measured 
data, this slow evolution is removed from the experimental data by applying a temporal 
smoothing of the signal, which cuts off all frequencies below 800 Hz. For consistency, the 
same smoothing is applied to BOUT data. 

A. Fluctuations: temporal and spatial characteristics, PDF 

The analysis of BOUT results shows that the h^at fluctuation amplitudes in LAPD data 
and the simulations have similar radial location near the cathode edge (r ~ 28 cm), where 
the background density gradient is largest, as shown in Fig. |6j The absolute values of the 
fluctuations are of the same order of magnitude, with simulated amplitudes smaller by a 
factor <2 than the experimental data. 

The comparison of the frequency power spectrum of the density fluctuations 5n/n to the 
LAPD measured spectrum is presented in Fig. [7J The spectra are integrated over the volume 
0.22 < r < 0.28 m, using a sliding Hanning window for averaging between the different 
simulation runs. Note that the total power in each spectrum is normalized to obtain the 
best fit to the experimental data. The power spectral shape from the BOUT simulation 
(for a range of v in values) is in a relatively good agreement with the measured spectrum 
(Fig. [7J b,d,e). At higher frequencies, > 10 kHz, the simulated spectra fall off faster than the 
measured spectrum. More studies are required to analyse the effect of additional features of 
the physical model (temperature profile and perturbations, sheath boundary conditions, etc.) 
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FIG. 6. (Color Online) Radial distribution of the average ion saturation current fluctuations 
normalized to the maximum equilibrium profile value. 

on the power spectrum, as well as the numerical resolution. High frequencies corresponding 
to smaller spatial structures are potentially more affected by finite resolution effects. As a 
check of numerical convergence in terms of box size, the spectrum for v in /Vt ci = 2 x 10~ 3 is 
calculated with four times the size of the azimuthal extent of the grid and the number of 
azimuthal grid points (N z = 128). The spectrum shape is similar to the original calculation, 
with a smaller grid size (Fig. [7J b and c). These simulations were performed without the 
explicit ion viscosity term in the vorticity equation ([3]), therefore the only viscosity is due 
to the numerical discretization. Inclusion of the ion-ion viscosity term as discussed in the 
Appendix |B] corresponding to the ion temperature Tj = O.leV does not significantly change 
the spectrum shape. The effects of viscosity in BOUT simulations is a subject of ongoing 
work. It is interesting to note that both the experimental and the simulation spectra exhibit 
an exponential power spectrum at higher frequencies (straight line on the log-lin plot), which 
is consistent with the presence of coherent structures^. 

Another important characteristic of the turbulence is the probability distribution function 
(PDF) of fluctuation amplitudes. The PDF of 5n/RMS(5n) fluctuations from LAPD probe 
data is integrated over a volume of plasma 0.22 < r < 0.28 m and compared with the 
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FIG. 7. (Color Online) Frequency power spectrum of the density fluctuations: LAPD measurements 
(a) and BOUT simulations (b-f) for b) u in /Q. ci = 2 x 10" 3 , c) v in /Q. ci = 2 x 1CT 3 with N z = 128 
azimuthal grid size, d) Vi n /Q, c i = 1 x 1CU 3 , e) Vi n /Sl c i = 2 x 10~ 4 , f) Vi n /Sl c i = 2 x 10~ 3 with ion 
viscosity at Tj = 0.1 eV. Experimental density profile, T e = 5eV. 

PDF from BOUT simulations in the same volume (Fig. [8]). There are no normalizations 
or fit factors involved in this comparison. The experimental and the simulated PDFs are 
similar, with the average relative fluctuation (\5n/n\) of 0.16 for the measured data and 
0.09, 0.09, 0.08 for BOUT simulations with v in /VL ci = 2 x 10~ 4 , 1 x HT 3 , 2 x 10~ 3 . The PDF 
for Uin/Qd = 2 x 10~ 3 case is the closest to the experimental data, which is consistent with 
the estimate of the neutral density in LAPD. Note that the distribution is mostly symmetric 
here because it is integrated over a radial inteval where the skewness is relatively low (Fig. [9]). 

Intermittent turbulence is observed in the edge plasmas of many experimental devices. 
This intermittency is usually attributed to generation and transport of coherent filaments of 
plasma, "blobs" or "holes"—. One signature of the presence of these structures is the non- 
zero skewness of the density fluctuation PDF. Typically, positive skewness associated with 
convective transport of "blobs" is observed in LAPD measurements in the region outside 
of the cathode edge (r > 28 cm). Smaller negative values, associated with the "holes" are 
usually observed inside the cathode radius. The radial profile of the skewness of Sn is shown 
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FIG. 8. (Color Online) Probability distribution function of 5n/RMS{5n) fluctuation amplitude in 
BOUT simulations and LAPD data. PDF is volume averaged in the interval 0.22 < r < 0.28 m. 



in Fig. O Except for the edge of the simulation domain which is affected by the imposed 
boundary conditions, the trend of the skewness profile, as well as the absolute values, is 
similar in BOUT simulations and in the LAPD data. 



Two-dimensional turbulent correlation functions are measured in LAPD using two probes: 
a fixed reference probe and a second probe that is moved shot-to-shot to many (~ 1000) 
spatial locations in a 2D plane perpendicular to the magnetic field. The reference probe 
remains at a fixed position that is close enough to the moving probe in the axial direction 
so that the parallel variation of the turbulent structures can be neglected. This allows to 
obtain the 2D spatial correlation function in the azimuthal plane. A similar "synthetic" 
diagnostic to post-process BOUT simulation results is constructed by calculating the cor- 
relations between a reference location and all other points in each azimuthal plane. The 
correlation length in BOUT simulation is of the same order, but larger than the measured 
value (Fig. HO). 
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FIG. 9. (Color Online) Skewness of 5n distribution in BOUT simulations and LAPD measurements. 
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FIG. 10. (Color Online) Correlation function for Ii ySa t fluctuations measured using a moving probe. 

B. Fluxes and sources vs. inferred source/sink in LAPD 

An inferred particle source is required to maintain the density profile close to the exper- 
imental values, as described in section IIIIC1 In BOUT simulations, the calculated source 
function S(r) that produces steady-state turbulence with the desired density profile ap- 
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pears to be positive inside r ~ 28 cm, and negative outside (Fig. (ITT]) ). Remarkably, the 
qualitative form and magnitude of the S(r) profile is consistent with the assumption that 
within ~ r there is an ionization source, and outside of ~ r there is a sink due to paral- 
lel streaming to the end walls. In LAPD, the field lines inside r ~ 28 cm connect to the 
cathode that produces the primary ionizing electron beam. The magnitude of the inferred 
source is close to the estimated ionization source and parallel losses sink for LAPD plasma, 
Ssource ~ n n n e < ov >j~ 2 x 10 20 m~ 3 /s and S S mk ~ niC s /L\\ ~5x 10 20 m~ 3 /s. 
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FIG. 11. (Color Online) Inferred particle source required to maintain the measured density profile. 



V. DISCUSSION 

A. Plasma transport 

One can calculate the effective diffusion coefficient by dividing the radial flux by the 
gradient of the equilibrium density, 

r 



D, 



eff 



VA, 



(11) 



,n 



The turbulence- driven radial particle flux in BOUT simulations peaks near the maximum 
density profile gradient and is close to diffusive model value, with the effective diffusion 
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coefficient -D e // ~ 3 m 2 /s on the order of Bohm value, DBohm, ~ 8 m 2 /s. The profiles of the 
average radial flux measured in LAPD and calculated in BOUT are shown in Fig. ffT2|) . The 
simulated flux is significantly lower than the measured flux in LAPD, but this difference 
is consistent with the difference in fluctuation amplitude shown in Fig. El It should also 
be noted that the measured flux profile is valid only for r > 28cm as the measurement of 
azimuthal electric field fluctuations in LAPD is affected by fast electrons generated by the 
plasma source^. 

While comparisons are made here to a diffusive model and particle transport in LAPD 
has been shown to be well modeled by Bohm diffusion^, the simulated plasma turbulence 
has features that are usually associated with intermittency, as discussed in section IIVA1 
non-Gaussian fluctuations, as seen by non-zero standard statistical moments. The role of 
coherent structures and convective transport in BOUT simulations is the subject of ongoing 
investigation. 
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FIG. 12. (Color Online) Left: Radial profiles of average n, and drii/dr. Right: radial particle flux 
r from the measured LAPD data and BOUT simulations. 



B. Numerical diagnostics and (E r ) equation 

An important part of the simulation effort is the framework for various diagnostics of 
the numerical solution. An obvious test of the solution is the check of the conservation 
laws (particle number, momentum, etc). However, depending on the choice of the radial 
boundary conditions and the form of the sources, the total number of particles, for example, 
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is not necessarily conserved during the simulation. A more appropriate diagnostics in this 
case is the local check of the balance of the terms in each of the time evolution equations 
(P~H2D, at each point in space and time. As well as ensuring that the solution is correct, 
calculation of this balance also provides a useful insight into the relative importance of the 
physical terms. 

This test can be done directly, using the explicit form of the equations ([T]|3]) as they appear 
in the code, or indirectly, by calculating the conservation of physical quantities not directly 
solved for by BOUT. For the direct tests, the balance of the terms in BOUT equations is 
satisfied to the machine accuracy when the appropriate finite difference schemes are used in 
the diagnostic module. 

Tests involving equations not directly solved by BOUT can be more subtle. An example 
is the equation for the azimuthal momentum. BOUT simulations in LAPD geometry show 
generation of self-consistent axisymmetric component of the potential. The dynamics of this 
zonal flow component can be derived from the vorticity equation as shown in Appendix [A] 

d t (NVg) = K) + ~ /Nd e ^-\ - V m (NVg) - fid r (w) (12) 

This expression is similar to the equation for the zonal flow component of the radial elec- 
tric field that can be derived directly from the azimuthal projection of the ion momentum 



equation 44 : 



^d r (r 2 (V r Vg)) - v m (V e ) + y, (d 2 T (V e ) + \ (V 6 ) - ^ (Vo)) 



d * W = ~^ d r ( r2 (VrVe)) - Vin (V e ) + n [d 2 r (V e ) + -d r (Vg) - - 2 (V e ) j (13) 

The first term in the right-hand side of Eq. f fl2|) and Eq. f fl3|) is the turbulent Reynolds stress. 
If a stationary turbulent state exists, the time average of the Reynolds stress, which is the 
driving term for the zonal flows, is balanced by the ion-neutral collisions and the ion viscosity 
terms. This equation is not derived directly from the vorticity equation, therefore BOUT 
solution satisfies it only to the extent that the underlying assumptions in the derivation of 
the vorticity equation ([3]) are fulfilled. 

Applying this diagnostic to BOUT output, it is found that the balance is well satisfied 
(within ~ 5%) for Eq. ( 1T2|) (with a small correction due to the linearization used in the 
inversion of Eq. (jlj)) and Eq. ( II 3p with a particular choice of the finite difference scheme 
in the advection operators - the fourth order central scheme. However, for the first order 
upwind scheme used in most of nonlinear simulations presented here, the balance is not 
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sufficiently well satisfied, which indicates that a higher numerical resolution is required to 
reach convergence for this diagnostic measure. 

At present, a quantitative match between the average (E r ) in the code and in the ex- 
perimental data has not been obtained. However, at present the model does not include 
all physics (e.g., sheath, biasing, temperature perturbations) that is certainly important for 
setting the average radial electric field. Improving the model by adding to it the missing 
physics to address matching of (E r ) is the subject of ongoing research. 

C. Future work 

A detailed verification study of linear instabilities in LAPD using BOUT 33 combined 
with a good qualitative and even partially quantivative agreement between nonlinear BOUT 
simulations and LAPD measuments presented here provides confidence in the relevance of 
these simulations to LAPD. It is, however, not a fully consistent first-principle model at 
present, and some potentially important physics is yet to be included. Most importantly, 
the experimentally measured temperature and flow profiles need to be matched and an 
evolution equation for temperature fluctuations, already implemented in BOUT, could be 
employed. The addition of a temperature gradient to the simulation would likely increase 
instability drive and result in a larger saturated turbulent amplitude and particle flux, more 
consistent with observation. Another significant improvement of the model that remains to 
be implemented is the sheath boundary condition at the end plates in the parallel direction. 
The axial boundary conditions are expected to be important in the formation of the radial 
electric field profile in LAPD, and a more physical description might help to understand and 
model the dynamics of the azimuthal flows and experimentally relevant potential profiles. 
The role of the ion viscosity on these flows is also a subject of ongoing work and has to be 
investigated in more detail. Additionally, including sheath boundary conditions along with 
temperature gradients, electron temperature fluctuations and an equation for temperature 
evolution can give rise to modifications to drift instabilities and introduce new modes, such 
as entropy waves 4 ^^ and conducting wall modes 47 . 

Although the linear calculations indicate that electromagnetic effects do not significantly 
affect drift wave instability in LAPD, magnetic field perturbations are required for Alfven 
wave studies. Alfven waves represent an important part of LAPD research, and nonlinear 
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simulations of Alfven waves using BOUT can contribute to the understanding of electro- 
magnetic turbulence in LAPD. 

VI. CONCLUSIONS 

A numerical 3-D model of plasma turbulence is applied to LAPD. The physics model 
includes equations for plasma density, electron parallel momentum, and current continuity, 
for partially ionized plasma. The model is implemented in the numerical code BOUT that 
is adapted for cylindrical geometry. This model has previously been successfully verified for 
a range of linear instabilities in LAPD 33 . 

Two different methods for average profile control are applied in the simulations presented 
in this study. One approach consists of suppressing the azimuthal average of the density 
fluctuations by directly subtracting it from the time evolution equation. The second method 
uses a source/sink radial function that is constructed to balance steady-state radial trans- 
port flows. Both methods successfully maintain the average density profile close to the 
experimental value, required for comparisons with the measured turbulence characteristics. 

Nonlinear BOUT simulations demonstrate a self-consistent evolution of turbulence and 
self-generated electric field and zonal flows, that saturates and results in a steady turbulent 
state. The simulated fluctuation amplitudes in the steady state are within a factor of 2 of 
the measured data, with a similar radial location near the cathode edge. The probability 
distribution function of the fluctuation amplitudes is comparable to the experimental dis- 
tribution. Statistical properties of edge turbulence, such as the skewness, are often used as 
indication of the turbulence intermittency 7 . In tokamak edge, the skewness is positive in 
the SOL (i.e. dominated by large amplitude events), but is sometimes negative inside the 
separatrix or limiter radius (i.e. with density "holes"), e.g. seen in DIII-D^ and NSTX— . 
This feature is similar to LAPD data, and reproduced in BOUT simulations. Despite the 
intermittent character of turbulence, as indicated by non-Gaussian PDF, the turbulent par- 
ticle flux magnitude is consistent with diffusive model with diffusion coefficient on the order 
of Dsohm- The inferred particle source/sink function required to maintain the simulated 
steady-state density profile close to experimental value is consistent with the estimates of 
ionization sources and parallel losses in LAPD discharge. 

The spatial and temporal structures of the fluctuations are consistent with the LAPD 
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measurements, however the correlation length in BOUT simulations is larger than in the 
experiment. 

Although some elements of the physical model still remain to be implemented (sheath 
parallel boundary conditions, azimuthal flow matching, electron temperature fluctuations, 
etc.), the agreement between certain features of the experimental data and simulations based 
on this relatively simple model lends confidence in the applicability of these simulations to 
the LAPD plasmas. 
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Appendix A: Azimuthal momentum equation 

The equation describing the evolution of the surface-averaged azimuthal momentum can 
be derived from the vorticity equation: 

d t zu = -v E ■ Vw - V||(iVu||) + b x ViV • Vv £ 2 / 2 - u in w + imV\w (Al) 

where the potential vorticity is 

w = V ± • (NV \<j>) (A2) 

Define surface and volume averaging by 

</(r ' M)) = ^Fi/ V fdzde (A3) 

</(r, 0, z)) v = ± f (f(r>, 9, z)) r'dr' (A4) 



Note the identity 



Vl/> = Vi (/) (A5) 
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For LAPD geometry, the convention is B z = —B p , so bo = — z. In normalized variables, 

v E = b x V0: 

v = _ hz \®± = \d± 

2 r 80 r 80 



or or 



(A6) 



The equation for the evolution of the azimuthal flows can be obtained from the vorticity 
equation ( 1A1I) by volume averaging. Applying Gauss theorem and assuming that the bound- 
ary conditions on the internal boundary are such that all surface integrals over the internal 
surface r = r a vanish, the following expression for the volume- average of the vorticity is 
obtained: 



w) v = — (Nd r 4>) = -— (NV e ) (A7) 



The ion-ion viscosity term: 



/„o \ 2ttL .„ . 2ttL „ . . . , . 

<Vi^) y = -y (d r w) = —d r (w) (A8) 

The advection term can be rewritten as: 

(v B • Vw) v = (V • (ot £ ) - tuV • v £ ) y (A9) 

where in straight field the last term vanishes. Applying Gauss theorem, the full divergence 
becomes a surface average: 

(V • (ujv e )) v = y (wV r ) (A10) 

The fourth term in Eq. (1A1I) can be written as a full divergence: 
bxViV- V— = V • ( — b x ViV J - — V • (b x ViV) = V • ( ^^b x VN J (All) 
The volume integral is then transformed info surface average 

bxViV.V^\ =M/Zi^ bxV /V.f\ (A12) 



V r\ 2 I V r\ 2 

Collecting all terms, one obtains the surface- averaged azimuthal momentum evolution 
equation: 

dt (NV e ) = K) + J /Ndg^-\ - v in (NV e ) - fid r (w) (A13) 
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Appendix B: Ion viscosity 

Including the ion-ion viscosity effects results in an additional term in the vorticity equa- 
tion 03]) • Viscous force V • II (where II is the stress tensor vv — t> 2 /31) in the ion motion 
equation induces a perpendicular drift v^ = b x (V • II) to the lowest order, which trans- 
lates into an additional perpendicular current component jj_ „ = enb x (V • II). To simplify 
the final form of the viscous term, it is assumed that density perturbations are small and the 
equilibrium gradients are negligible compared to the perturbed quantities. Then the extra 
term due to viscosity in the vorticity equation can be written as 

V • j Ltll = enV ■ {b x (V • n)} = -enb ■ {V x (V • II)} (Bl) 

Using Braginskii expressions for the stress tensor 36 with perpendicular ion velocity v i± = 
\e = b x V0 to the lowest order, it can be shown that the extra term in the vorticity 
equation ([3]) is ii^S\w. 

Depending on the ion magnetization, one should choose either the magnetized or unmag- 
netized viscosity expression for the coefficient \io \i^ — r]\ — 0.3nTj/f^Tj for Vt ci Ti > 1 or 
Hi — tJq — 0.96nTiTi for f2 C jTj <C 1. The estimate of the ion magnetization parameter f2 C jTj 
for typical LAPD values (He 4 , B = 400 G, m ~ 3 x 10 18 m" 3 , T t ~ 1 eV) is close to unity. 
The ion temperature in LAPD is not directly measured; the estimate from the electron-ion 
energy exchange and the parallel losses balance indicates that T, is in the range 0.1 — 1 eV. 
Preliminary BOUT simulations with ion viscosity are consistent with the experimental mea- 
surements for low ion temperatures, Tj ~ 0.1 eV, which corresponds to unmagnetized ion 
regime. 

Appendix C: Discretization in the parallel coordinate 

Consider the local drift mode dispersion relation for the simplest case, without electron 
inertia and electromagnetic terms, when it becomes a quadratic equation, as given in ele- 
mentary plasma textbooks 50 : 

(w- l)zcj||+u; 2 = 0, (CI) 



where 



rv|| \ ii c jii ce 



k±J Vei^* 
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u* = k±v pe = -±— g- (C3) 

and uj is normalized to w*. 

Now consider the effect of finite-difference discretization on Eq. (ICljl . For simplicity 
assume no radial structure so that the radial part of the solution can be dropped. There are 
just two coordinates then - the drift wave propagation direction y and the parallel direction 
z. 

Now let's focus on the parallel discretization. Parallel derivatives that are represented 
by iku in the Fourier form will become something different in the discretized equation, 
depending on the type of discretization. For example, applying the 2 nd central difference for 
a single mode exp(ik\\z) one obtains 

df exp(ik\\Zj + i) — exp(ik\\Zj_i) ism(k\\h) 

Tz "* 2h = h (C4) 

Here h is parallel grid spacing, Zj=jh. 

One can notice that at large wavenumbers, k ~ 7t/h, the finite-difference representation 
is very poor. In this example, from Eq. IC14 for a\\ ^>1 the growth rate scales as 7 oc 1/fcn, 
i.e., large k\\ should stabilize the modes. Conversely, in the discretized dispersion relation 
it will become 7 oc l/(sin(k\\h)) 2 which would become singular at the Nyquist wavenumber 
k\\ = 7t/h, which can be manifested in unphysical behavior of such modes. The possibility 
of unphysical behavior of high-k modes caused by spatial discretization, in particular the 
"red-black" numerical instability, is well-known in the CFD community, and historically the 
main remedy was using staggered grid^i. A more recent popular method is discretization on 
collocated grids adding a dissipative biharmonic (i.e. 4 th derivative) term to suppress the 
"red-black" instability (e.g., the Rhie-Chow interpolation). 

Consider the effects of staggered grid for the discretized drift mode dispersion relation. 
Assume that Ni, 0, w are specified on one grid, and j||, V\\ e are specified on another grid 
shifted by h/2. Then, k?, in the dispersion relation (jCll) can be tracked down to the deriva- 
tives 

9j\\ ,. . w , / ., , ,-.s i sin(kh/2) ,_,_, 

~§ ->■ (l\\,J-J\\,i-i)/h = exp(-ikh/2) ±j^J. (C5) 



and 



d(£>\\ , , ,, , , , , , isin(kh/2) ,„ . 

_m ^ ( 0i+1 _ f.yh = exp ( ik h/2) ^L2, (C6) 
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which combine to produce 

One can note that Eq. (1C7|) does not become zero for any mode supported by the grid, 
—7i/h < k < n/h, which is an important improvement. 

Since the staggered grid approach is more cumbersome for implementation, in particular 
for parallel computation, it is desirable to stay with collocated grids, if possible. In this 
example one can come up with discretization Eq. (1C5HC6J) on collocated grids by combin- 



ing right-sided and left-sided I s * order discretization. This method, which we call "quasi- 
staggered grid" , has been successfully applied in BOUT to eliminate spurious modes due to 
parallel discretization. 
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